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Abstract 

We demonstrate an efficient nonequilibrium Green's function transport calculation procedure 
based on the real-space finite-difference method. The direct inversion of matrices for obtaining the 
self-energy terms of electrodes is computationally demanding in the real-space method because the 
matrix dimension corresponds to the number of grid points in the unit cell of electrodes, which is 
much larger than that of sites in the tight-binding approach. The procedure using the ratio matrices 
of the overbridging boundary-matching technique [Phys. Rev. B 67, 195315 (2003)], which is 
related to the wave functions of a couple of grid planes in the matching regions, greatly reduces 
the computational effort to calculate self-energy terms without losing mathematical strictness. In 
addition, the present procedure saves computational time to obtain Green's function of the semi- 
infinite system required in the Landauer-Biittiker formula. Moreover, the compact expression to 
relate Green's functions and scattering wave functions, which provide a real-space picture of the 
scattering process, is introduced. An example of the calculated results is given for the transport 
property of the BN ring connected to (9,0) carbon nanotubes. The wave function matching at the 
interface reveals that the rotational symmetry of wave functions with respect to the tube axis plays 
an important role in electron transport. Since the states coming from and going to electrodes show 
threefold rotational symmetry, the states in the vicinity of the Fermi level, whose wave function 
exhibits fivefold symmetry, do not contribute to the electron transport through the BN ring. 

PACS numbers: 71.15.-m, 72.10.-d, 72.80.Rj, 73.40.-c 
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9 I. INTRODUCTION 



10 Recently, electronic-structure calculations have become an important tool for investigat- 

11 ing the physics and chemistry of nanoscale systems with the miniaturization of electronic 

12 devices. The electron-transport properties of nanoscale systems have been studied actively 

13 because they are of significant importance from both fundamental and practical points of 

14 view. Owing to the complexity of the problem, such studies are strongly dependent on the 

15 existence of reliable numerical treatments based on first-principles approaches. 

16 A number of first-principles calculation methods for the electron-transport properties of 

17 nanoscale systems have been proposed so far. They are roughly categorized into two ap- 
is proaches. One approach uses the nonequilibrium Green's function (NEGF). The relation be- 

19 tween conductance and Green's function has been derived within the nonequilibrium Keldysh 

20 formalism.- This approach has been used extensively in connection with tight-binding mod- 

21 els and first-principles methods employing localized basis sets consisting of either atomic 

22 orbitals^"- or Gaussians.^ In this formalism, by making use of Green's functions with en- 

23 ergies of nonreal numbers, the electronic structures of the isolated states in the transition 

24 region, which are not easily treated only by the energies of real numbers, can be included into 

25 the total charge density of the system. On the other hand, the NEGF method has not been 

26 employed with the real-space^"- and plane-wave formalisms so far because the large num- 

27 bers of grid points or plane waves prevent the direct inversion of N x x N y x iV z - dimensional 

28 matrices to obtain the surface Green's functions and self-energy terms of electrodes, where 

29 N x , N y , and N z are the numbers of grid points along the x-, y-, and ^-directions in the unit 

30 cell of the left or right electrode and electrons flow along the ^-direction [see, e.g., Fig. [Tj. 

31 The other approach is to compute the scattering wave functions from which transmission 

32 coefficients can be obtained. In addition, the scattering wave functions provide a direct 

33 real-space picture of the scattering process. This approach has been employed by combining 

34 it with techniques where real-space grids and/or plane wave basis sets^"— are used to de- 

35 scribe wave functions and potentials. The easiest way to obtain scattering wave functions is 

36 to solve the Lippman-Schwinger equation where the semi-infinite electrode is replaced with 

37 a uniformly distributed charge background, i.e., "jellium."- The scattering wave function 

38 can alternatively be calculated by the wave-function-matching approach proposed by Fuji- 

39 moto and Hirose, i.e., the overbridging boundary- matching (OBM) method.- 1 ^ In the OBM 
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40 method, several parts of Green's functions of the transition region are computed to set up 

41 the wave-function- matching formula and there is no limitation of the usage of jellium, i.e., 

42 more realistic atomic and electronic structures of the electrode can be taken into account. 

43 Moreover, the effect of electrodes is included in the matching formula as N x x N y x Mf- 

44 dimensional ratio matrices, where A/} is the order of the finite-difference approximation^ for 

45 the kinetic energy operator in the Kohn-Sham equation and is much smaller than N z . 



46 In this paper, we propose the real-space finite- difference (RSFD) NEGF scheme and 

47 demonstrate that the surface Green's functions and self-energy terms of the electrodes re- 

48 quired in the NEGF method can be obtained from the ratio matrices in the OBM method. 

49 The dimension of the matrices for the surface Green's functions and self-energy terms that 
so are directly inverted is reduced from N x xN y x N z to N x x N y xMf in the present scheme while 

51 keeping the rigorousness of the mathematical formulation. This advantage also saves com- 

52 putational effort to calculate Green's functions in the transition region because the number 

53 of elements of the Green's-function matrix required for the transport calculation is propor- 

54 tional to the dimension of the matrices of the self-energy terms. In addition, we prove that 

55 scattering wave functions, which help us to interpret the scattering process and are usually 

56 calculated by wave function matching methods, can be obtained within the framework of 

57 the RSFD NEGF method. 



58 The rest of this paper is organized as follows: in Sec. II, we relate the important quan- 

59 tities of the NEGF method to those of the OBM scheme, and introduce the RSFD NEGF 

60 method. In Sec. Ill, we present an example showing that scattering wave functions help 
ei us to interpret the scattering process by computing the transport properties of a BN ring 

62 sandwiched between carbon nanotube electrodes (C/BNNT). In Sec. IV, we summarize our 

63 procedures. Finally, a mathematical proof is given in appendix. 

3 



64 II. NEGF METHOD USING THE COMPUTATIONAL PROCEDURE OF THE 

65 OBM METHOD 



66 A. Green's function of a whole system including the transition region and two 

67 semi-infinite electrodes 

68 Let us consider Green's function of a system composed of the transition region sandwiched 

69 between two semi-infinite crystalline electrodes, as shown in Fig. [TJ within the framework of 

70 the RSFD scheme. Two-dimensional periodicity in the x- and ^/-directions is assumed and a 

71 generalized z-coordinate Q instead of z\ is used, because a couple of grid planes are involved 

72 in wave-function and Green's-function matching when higher-order finite-difference approx- 

73 imation is employed (see Fig. [TJ). The exchange- correlation effect is treated by the local 

74 density approximation^ or generalized gradient approximation 1 ^ of the density functional 

75 theory.—^ 

76 The Green's-function matrix involves the inversion of an infinite matrix corresponding to 

77 the Hamiltonian matrix of the whole system H(ku), where k\\ is the lateral Bloch vector. 

78 As shown in Fig. [3J we are, however, interested in the finite part of the Green's-function 

79 matrix, 







Blt 





#(fc„) = 


Sir 


H T (k\\) 


Btr 







B ] 

TR 





so where the borders of the partitioning of H(ku) are drawn according to the dashed lines in 
si Fig. [TJ the submatrix Hj-^ku) contains the matrix elements in the transition region, HT(k\\) 

82 (Hu(k\\)) corresponds to the semi-infinite left (right)-electrode region, and Blt (Btr), which 

83 is nonzero only for some connection area between the transition region and the left (right) 

84 electrode, is the interaction between the transition region and the electrodes. Hx(k\\) in the 

85 transition region is treated as a general nonsparse matrix here and in subsequent subsections, 
se On the other hand, Hi,{k\\) (Hn[k\\j) is a block tridiagonal matrix in the RSFD formalism, all 
87 of the block-matrix elements of which are N(= N x x N y x A//) dimensional. In addition, B^t 
as (Btr) is a zero matrix except for one iV-dimensional block-matrix element B((_i) (B(( m+ i)), 

89 as illustrated in Fig. |2j In practice, Aff corresponds to the number of x-y grid planes involved 

90 in the function- matching region Q since the order of the finite-difference approximation is 
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chosen so as to include the nonlocal region of pseudopotentials in the matching region^ as 
well as to obtain sufficiently accurate results with the used grid spacing. 
Green's function of the whole system is defined as 



G(Z,k\ 



H(k\ 



G L (Z,k l{ ) 


G LT {Z,k\\) 


G LR (Z,k\\) 




G T (Z,k\\) 


G TR (Z,k\\) 


G RL {Z,k\\) 


G RT {Z,k\\) 


G R {Z,k\\) 



94 where Z{= E + irj) is a complex energy variable. From the matrix equation 
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95 that is, 



G LT (Z,k\\)G T (Z,k\ 



B\ T G LT {Z,k\\)+ Z-H T {k\\) G T (Z,k 



Z-H L {k\\) 
■ B TR G RT (Z, k 

Z-H R (ku 



Blt 
)=I 



B, 



t 

TRi 



(2) 



(3) 



(4) 



G RT (Z,k\\)G T {Z,k\\)-' 

one sees that Green's function of the whole system Gx{Z^k\\) can be portioned to the 
transition region as 



G T (Z, fen) =[Z- # T (fc||) - J2 L (Z, fc||) - £ B (Z, fc||)j . (5) 

98 Note that Eq. (JSJ) is equivalent to Dyson's equation in the standard form^S, of 

g t (z, fen) = g T (z, fen) + g T (z, fen) [j: l (z, fe„) + £ R (z, fc,,)] 6 T (z, fc,,). (6) 

99 Here, ^2 L (Z,k\\) and ^2 R (Z,k\\) are the self-energy terms of the left and right electrodes 
wo defined by 



± L {Z,k\\) = Bl T Q L {Z,k\\)B LT 
J2 R (Z,k l{ ) = B TR g R (Z,k {l )P TR , 



(7) 



ioi where 



G L {Z,k\\ 

Q R {z,k\ 



Z-H L {k\\) 
Z-H R {k\\) 



(8) 



102 are Green's functions of the semi-infinite left and right electrodes with right- and left-side 

103 truncations, respectively. In addition, Qt{Z, ku) is the Green's function associated with the 

104 isolated transition- region Hamiltonian Hx(k\\): 



-1 



Z-H T (ku) . (9) 



105 We used the script capital letter Q for describing Green's function of a semi-infinite sys- 
ioe tern with one-side truncation as well as an isolated (two-side truncated) system to prevent 
107 confusing [Z — Ha] with Ga defined by Eq. (j2J), where A = L, T and R. From Eqs. §5§ 
we and (jBJ), one sees that Q T (Z,ku) is extended to Gt(Z, ku) so as to include the effects of 
109 semi-infinite electrodes through the self-energy terms Y1{lr}(Z, k\\). 



no B. Evaluation of self-energy terms 

hi This subsection is devoted to the evaluation of the surface Green's functions of the left- 
n2 and right-electrode regions G{l,r}(Z, k\\) and the self-energy terms J2^ L m(Z, ku). Hereafter, 
n3 we omit the branch of the lateral Bloch vector k\\ for simplicity To set up the Hamiltonian of 
H4 the electrodes, the Kohn-Sham effective potential is obtained using the unit cell consisting of 
us N x x N y x N z grid points under the periodic boundary condition [see Fig. [3]. The Green's- 
n6 function matrices for electrodes are computed using the recursive technique proposed by 
in Lopez Sancho et alr^ or Guinea et al.— Alternatively, the matrices are evaluated by solving 
us a quadratic eigenvalue problem.— In both schemes, the peculiar characteristic that the 
n9 N x x N y x ^-dimensional matrices of the Hamiltonian of the electrodes are the same for 

120 each unit cell of the electrodes is employed, which implies that the dimension of the matrices 

121 for Green's functions and the self-energy terms becomes N x x N y x N z . When the number of 

122 grid points increases, the computation of these matrices is demanding. In the RSFD NEGF 

123 scheme, since Blt {Btr) has only one nonzero N(= N x x N y x A/})-dimensional block- 

124 matrix element JB(£_i) (B(( m+ i)) [see Eq. (JTJ and Fig. |2j, the self-energy terms, which are 

125 calculated by Eq. (J7|), are found to take the very simple form of 
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El(S) 



El((o;Z) o •■• o 

•■■ 



■■■ 

o ••• o E R (Cm+i;Z) 



(10) 



where 



E L (Co; z) = B(uyg L (C-x, C-i; ^)£(C-i) 

X^Cm+i; ^) = ^(Cm+l)^ii(Cm+2, Cm+2! Z)B((, m+1 ) ] (H) 

with ^{L,ii}(Cit, 0; Z) being the iV- dimensional (A;, I) block-matrix element of G{l,r}(Z). In 
practical calculations, A// is much smaller than N z . 

To obtain Green's function of the whole system G(Z), it is sufficient to calculate N- 
dimensional matrices for the self-energy terms. However, although the Hamiltonians of the 
Mth and M + 1th unit cell are identical, the sliced matrices of the Hamiltonian of the 
electrodes, #(C/ M ) and 1 ), are not the same, where H(Q ) is the lih iV-dimensional 

diagonal block-matrix element of the N x xN y x iV z -dimensional Hamiltonian matrices of the 
electrodes. Thus, we cannot use the two procedures above mentioned.— ~— We introduce a 
computational procedure for obtaining N- dimensional matrices of the surface Green's func- 
tions and self-energy terms using ratio matrices in the OBM method/ 4 which are computed 
by solving the generalized eigenvalue problem for the periodic bulk model shown in Fig. [3J 



n x (z) 



MC!~\z) 
<MCf + U) 



X n (Z)U 2 (Z) 



®n((l M+ \Z) 



(12) 



where 



MZ) 
Mz) 











e(cf, Cf ; z)B(C!Y e(cf, C; z)b& 



-An 



(13) 
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139 and O(Cjfe) C/ M j Z) is the AT- dimensional (A;, I) block-matrix element of Green's function of 

140 the truncated part of the periodic Hamiltonian in the Mth unit cell. In addition, Eq. f fl2|) is 
hi the analytically continued equation of Eq. (19) in Ref. Q- Note that the generalized eigen- 

142 problem of Eq. (fT2"j) suffers numerical error owing to the extremely large and small absolute 

143 values of A(Z) in some cases, which prevent us from conducting an accurate computation of 

144 the eigenstates {$ n }- To avoid this numerical difficulty, we introduce the following ratios of 
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145 the generalized eigenstates, which are proposed in Ref. 

R p (Co;Z) = Q p (C-i;Z)Q p (Co;Z)- 1 

R%Cm+2\ Z) = Q*(Cm +2 ; Z)Qi(( m+1 ; Z)~\ (14) 

146 where 

Q p (0; Z) = [$?(0; Z), $f(0; Z), & N (( i; Z) 

Q\Q- Z) = [$f (0; Z), Z), & N (& Z)] . (15) 

H7 Since the eigenvalues A n (Z)'s for a nonreal Z are divided evenly into two groups with |A„| > 1 
us and |A„| < 1,— we set up the N- dimensional matrix Q p (d; Z) and Q q (d; Z), which gathers 
H9 the N eigenstates {$£(0; Z)} and {^((f, Z)}, n = 1, 2, N, for a nonreal Z with |A n | > 1 

150 and |A„| < 1, respectively, using the solutions of Eq. ( fl2|) . Then the accuracy of ratio 

151 matrices is improved using the following continued-fraction equations in a self-consistent 

152 manner [see Eq. (25) of Ref. 

R p ((? +l ;Z) = e(C,C;Z)B(C) 

+e(C, Cf; z)B(Cy [iF(Cf ; z)- 1 - e(c 1 M , Cf; z^O 1 ] 1 ©(Cf, C; z)B{C) 

+0(Cf, z)5(C) [R^cr 1 ; zy 1 - 0(C z)B(C)] 1 e(e, Cf; z)5(O t . 

(16) 

153 Now, we prove the following relation that gives a definite description of the surface Green's 

154 functions (or self-energy terms) in terms of the ratio matrices of the generalized eigenstates 

155 of Eq. ( fl~2l) . The surface Green's functions of the left and right electrodes are explicitly 

156 expressed as 

^(C_i,C-i;Z) = i? p (Co;Z) J B(C-i)- 1 

G R (Cm+2, Cm+2; Z) = R% m+2 ' Z)B(( m+1 )^\ (17) 



157 and the self-energy terms Eq. (TTTj) are given by 

Z L (Co;Z) = B(C- 1 yR p (Co;Z) 

EhCOh-i; Z) = B(( m+1 )R«(( m+2 ; Z). (18) 

158 Hereafter, we concentrate on proving the surface Green's functions and self-energy terms 

159 of the left electrode because those of the right electrode can be derived in a similar manner, 
wo The proof is derived using the results reported by Lee and Joannopoulos:— Green's function 
«>! G{l,r}{0, Cm', Z 7^ real) of a crystalline bulk is a decaying function (i.e., G{l,r} — > 0) as 

162 \l\ — > oo with m fixed (or as \m\ — > oo with / fixed). The eigenstates are ^-independent 

163 functions with respect to Q with a decreasing or increasing property such that 

^ n {Ci-sL,Z) = {\ n )- s ^ n {CuZ) (19) 

164 Here, s is an arbitrary positive integer and L is an integer associated with the length of 

165 periodicity in the z-direction. For example, L = N z /Mj when N z is a multiple of J\ff. 
lee Furthermore, {<& n (Q;Z)} satisfies 

- s/_ 2 ^(0_ 2 ; z) + A-i^(0-i; z) - b^KCi; z) = o, (20) 

167 where 

A i = Z-H(C l ) and Bi = B(d). (21) 

168 Note that Eq. f[2Hl) is the Kohn-Sham equation when a complex number Z is replaced with 

169 a real number E. 

170 In the left electrode, the surface Green's function Gl(CiX-i'i Z) is expressed in terms of 
in QP(Cf,Z): 

G L (( h C_i; Z) = Q p (Cv, Z)Q?(( ; Zy^B^y 1 (I = -1, -2, ...,). (22) 

172 In the following, the derivation of Eq. fT22|) is demonstrated. It is straightforward from the 
m definition that {Gl((i,(-i\ Z)} satisfies 
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174 that is, 



-£i 5 £L 5 ,-l + A_30_4,-l " £-40-3,-1 = o 
-Bl 4 0_4,_i + A_ 3 0-3,-l " £-30-2,-1 = 
-Si 3 ^_ 3) -1 + ^-20-2,-1 " £-20-1,-1 = 

175 and 

- Bl 2 G-2,-i + A_xG-\-\ = I, (25) 

ne where 

£,-i = &(Ci,C-i;^)- (26) 

177 From the facts that the iV-dimensional Z ^ real) decays deep inside the left 

178 electrode (/ — > — oo),— we see that {$£(Cz; Z)} exhibits a linear independence of the decaying 

179 sequences, where n = 1,2, ...,N, and Eq. (12^1) are the same sets of simultaneous linear 
lso equations as Eq. (120]) for I < —1. Therefore, Ql is expanded in terms of {$£} and expressed 
lsi as 

" N N N 

/m*£(G; z), E /^(0; z), .... J] /niV^(0; 

n=l n=l n=l 

Q p (0;^)F (/ = -!, -2,...), (27) 



&(C/,C-i;^) 



182 where {f nn '} is a set of unknown expansion coefficients forming an N- dimensional matrix F. 

183 For simplicity, the dependence of f nn i and F on Z and k\\ is ignored. By inserting Eq. ( |2~T]) 

184 for Z = —1 and —2 into Eq. ( 125]) and subsequently using Eq. (120]) for Z = 0, we obtain 

F = QV(^Z)-\B_ 1 )-\ (28) 

18 5 By substituting Eq. QggJ) into Eq. ([27]), Eq. ([22]) is obtained. 

we Next, we carry out the limiting procedure in Eq. ( 122]) to obtain the retarded Green's 

187 function. Since the pairing characteristic of A(£)'s implies that there are always equal 

188 numbers of reflected waves $^ and transmitted waves <£>^ ra for a given energy E and lateral 

189 Bloch wave vector ku, we obtain 

lim <S>l{Q; E + ir))=$ r n e f {Cf, E), 
lim Q p (( l ;E + irj) = Q re f(( l ;E), 

r;-i>0+ 

lim RP(C l ;E + irj) = B re f(Cf,E), (29) 
10 



190 where Q ref {(i]E) is an iV-dimensional matrix that gathers left-propagating Bloch waves 

191 (lA(-E')l = 1 and Re(A(£')) < 0) and right ward increasing evanescent waves (|A(-B)| > 1): 



Q"'(G; E) = [^T f (Ci; £), $ r 2 e/ (0; E), ^ e/ (0; e)\ , (so) 

192 and 

iT e/ (G; E) = W e/ (0; E)~\ (31) 

193 Note that the eigestates with \\{Z) \ = 1 are absent when r] ^ 0,— while those with | \{Z) \ = 1 

194 exist in the case of real-number energy.— 

195 Finally, the retarded Green's function is 

g r L (C-i, C-i; e) = iim g L (u, c_ i; e + ir,) 

= Q^{^ 1 -E)Q^{C, 0] E)-\B_ l )- 1 

= R re f(( ;Z)B((_ l )- 1 , (32) 

we and hence the retarded self-energy term is 

ElW^O = Km E { L,ii}(^ + (33) 

197 Several researchers have investigated the representation of the surface Green's functions 
we by generalized eigenstates.- 1 ^ 1 ^ Note that the above relations are particularly attractive be- 

199 cause they allow us to directly evaluate the retarded surface Green's functions and retarded 

200 self-energy terms at a purely real E without using the finite broadening (or smearing) param- 

201 eter r]. More interestingly, we can reduce the dimension of the matrices of Green's functions 

202 and the self-energy terms from N x x N y x N z to N while keeping mathematical rigorousness. 



203 C. Evaluation of whole Green's function in the transition region 

204 Green's function of the whole system portioned to the transition region, Gt{Z), is given by 

205 Eq. (J5J). In this subsection, we present the analytic expression of Eq. ([2]), which is equivalent 

206 to the exact solution of Dyson's equation Eq. (jH]). In the general case of the Hamiltonian 

207 matrix H T being not necessarily sparse, the simplest way to calculate Gt{Z) of Eq. §5§ 

208 might be to carry out direct matrix inversion. In practice, however, it is computationally 

209 difficult to perform inversion calculations with large matrices. In what follows, we show that 

11 



210 there exists an efficient approach to computing the non-Hermit Gt(E) on the basis of the 

211 OBM scheme. 

212 Let us consider the Zth column of Gt(Z), i.e., [Gt (Co, Ci] Z), Gt(Ci,Ci] Z), • Gr(Cm+i, Cl'i Z)Y 

213 (/ = 0, 1, 2, m + 1). From Eq. (jSJ), one sees that the Zth column satisfies 

G T (Co,Ci;Z) 
G T (Cx,Ci;Z) 



z-h t -El(z)-E r (z) 




















I 


















G T (CuCuZ) = I <r- the hh (34) 



Gr(Cm, Oi Z) 
Gr(Cm+iXh Z) 

214 by virtue of a simple form of the self-energy matrices, Eq. f lTU]) . Using Green's function of 

215 the truncated part of the Hamiltonian Qt(Z) defined in Eq. (Q, the whole Green's function 

216 in the transition region is given by 



Gt(CoXuZ) 

G T (CiXi;Z) 



Gt(Ci,0;Z) 



GT(Cm, Ch Z) 

Gr(Cm+i, Ci]Z) 



Qt(Z) 



J2 L (C ;Z)G T (Ca,Cr,Z) 




/ 





Sfl(Cm+i; Z)G T (Cm+l,Ch z) 



<- the Ith 



(35) 



217 Equation ( 1351) is a matching relation with regard to Green's function {Gt(0c) 0)}; which 

218 is an analogous one with regard to the wave function {^(z/.)} [see Eq. (1) of Ref. Jjj]. 

219 The surface Green's-function matching theory has been pioneered by Garcia-Moliner and 

220 Velasco.— From Eq. (134")) . we see that once Green's function of the truncated part of the 

221 Hamiltonian Qt(Z) = [Z — Ht] 1 is known, the elements of the whole Green's function 
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222 Gt((o, Ci] Z) an( i Gr(Cm+i, Ci', Z) are calculated using 



£t(Co,Co; z)J2l((o; z) — i GriCoXm+i; z)J2 R {Cm+i; z) 

QriCm+i, Co', Z)J2l(Co] Z) Qr(Cm+l, Cm+1, Z)J2 R (Cm+l', Z) — I 



X 



G T ((a,(i;Z) 
GriCm+i, 0',Z) 



Gt({ ,(i;Z) 
Qt (Cm+i , 0; z) 



(36) 



223 Here, Gr(Ck, Ci', Z) is the iV-dimensional (k, I) block-matrix element of Qt(Z). Equation 

224 which is a 2N simultaneous linear equation with respect to Gt(Co, Ci', Z) and Gr(Cm+i, Cu Z), 

225 manifests boundary- value (surface) matching for Green's function of the whole system. To 

226 calculate the electronic structure in the transition region, the diagonal elements of the 

227 Green's-function matrix, Gt{Ci,Ci', Z), are required. On the other hand, when we are only 

228 interested in the transport property, it is sufficient to compute Green's functions on the 

229 matching planes, G r T (Co, Co; E), G T (C m +i, Co; E), G r T (Co, Cm+i! E), and G T {C m +i, Cm+u E). In 

230 the following, we show the analytic expressions for the solution of Eq. (136]) in the cases of 

231 / = and m + 1, which are required in the calculation of the conductance in the NEGF 

232 formalism. 

233 For I = and I — m + 1, 



G T (Co, Co; z) = g T (Co, Co; z) [i - £ L (Co; z)Q T (Co, Co; z) 
\ i _1 

Gr(Cm+i, Co; Z) = I - GriCm+l, Cm+1) Z)YjR(Cm+l, Z) 

x g T (C m+ i, Co; z) [i - £ L (Co; z)Q T (Co, Co; z) 
Gt{Co, Cm+i; z) — I - Q T {Co, Co; z)Y jL {Co', z) 

X I - J2 R (Cm+l', Z)Q T (Cm+l, Cm+1, Z) 
GriCm+l, Cm+lj Z) = QxiCm+l, Cm+1) Z) I — ^2 R {Cm+l', Z)QT(Cm+l, Cm+1, Z) 



Qt\Co, Cm+1, Z) 
-1 



i -1 



(37) 



234 where Qt is a modified Qt under the influence of the self-energy term X]{l,r}' which is 
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235 expressed as 



Gt{(o, Co; Z) = Gt(Co, Co', Z) 

+Gt(Co, Cm+i; Z)J2n(Cm+l] z) 



X 



I — GriCm+i, Cm+i; Zj^niCm+i] Z) GriCm+i, Co! Z), 

GriCm+l, Cm+li Z) = GriCm+l, Cm+li ^) 

+<MCm+i, Co; ^)Z1l(Co; z) 



X 



-f - £t(Co, Co; Z)J2 l (( ; Z) Qt (Co, Cm+i; Z). 



-1 



(38) 



236 It is easy to ensure that the Gts given by Eqs. (|37|) - (|38|) satisfy Eq. (|36|) . and therefore, 

237 they are the exact analytic solutions of Eq. (jSJ) as well as Dyson's equation of Eq. ([6]). 

238 Finally, a retarded Green's function is obtainable by carrying out the limiting procedure: 



G r T (Ck,Ci;E)= lim G T (C kl Ci;E + i v ). 

r/^0+ 



(39) 



239 Note that the block matrices of Green's function, G^(Co> Co; E), G^(C m +i, Co; E), G^(Coj Cm+i 

240 and G r T (C m +i-, Cm+ij E), are not A^. x N y x iV 2 -dimensional, but iV-dimensional because the 

241 matrices of the self-energy terms have already been reduced to N- dimensional in the pre- 

242 ceding subsection. 



243 D. Description of scattering wave function in terms of whole Green's function 

244 We next show that the relationship between the retarded Green's functions and the 

245 scattering wave functions is expressed as 



^•(0; E) = iG r T (Ci, Co; £)r L (Co; £)$f (Co; E ) 



(0 < / < m + 1) 



(40) 



246 where r^Cz; E) is the coupling matrix, which describes the 'coupling strength' of the tran- 

247 sition region to the left electrode at Co, and is defined by 



izCCo; E) = i EKCo; E) - EKCo; E) (E1(Co; e) = £2(Co; e?) . (41) 



248 From Eq. (6) of Ref. 



26l . the scattering wave function incoming from deep inside the left 
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249 electrode is expressed as 



%(Co;£) 



B(C-i)t*?(C-i; E) - El(Co; £)$r(Co; E) 










(42) 



250 The incident wave from the right electrode can be derived in a similar manner. By the 

251 definition of the retarded Green's function of the whole system in Eq. dSJ), Eq. (j4"2j) is 

252 rewritten as 



*i(Co;£) 
*i(Ci;^) 

*,-(Cm;^) 
^•(Cm+i;^) 



Gt(E) 



fl(C-i W (C-i; - EKCo; £)<*>f(Co; 



o 
o 



(43) 



253 Now the ratio matrix i? m in the left electrode (/ < 0) is introduced along a similar line 

254 into the definition of R re ^: 



255 where 



R m (( l ;E) = Q m (( l „ 1 ;E)Q m (( l ;E)-\ 



(44) 



(45) 



256 which is assumed to include not only ordinary right-propagating incident Bloch waves but 

257 also leftward- decreasing evanescent waves. From the definition of R m ((o,E), it is straight- 

258 forward to state that 

$f (C_i; E) = /T (Co; (Co; E). (46) 
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259 Furthermore, the relationship between the retarded self-energy term Y7l(Co] E)^ and i? m (Co; E) 

260 is expressed as 

J2 r L (Co;Ey = B(C-iYR m (Co;E) (47) 

261 similarly to that of Ez(Co5 E) and R ref {( ] E) in subsection [TEE From Eqs. flS}, (gSJ and 

262 (I47p . the scattering wave function ^(0; E) for < / < m + 1 can be written as 

^•(0; J5) = G^(0, Co; £0 [^(C-i)^f (C-i; e) - £ r L (C ; £)$f (Co; E) 
= -G r T (( h Co; £Q fc^Co; E) - EKCo; $f (Co; J5) 



= iG?.(Ci,Co;S)r L (Co;S)^(Co;£0. 

263 The derivation of Eq. ( 140"]) in a different manner is given in Re: 

264 N x x N y x ^-dimensional matrices is also introduced in Ref. 



(48) 

8|, and the expression using 
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265 E. Conductance 

266 We finally address the problem of electronic transport within the framework of the RSFD 

267 Green's-function approach. Here, we consider the case of the incident wave $* n incoming 

268 from deep inside the left electrode. We can prove that the Landauer-Biittiker formula^ 

269 G = 2e 2 /hJ2i j\tij\ 2v i/ v j describing the conductance G has the expression in terms of 

270 Green's functions G^f' a ^ and the self-energy matrices E{2'ii}> where e is the electron charge, 

271 h is Planck's constant, is the transmission coefficient for the jth incident wave to the 

272 ith outgoing wave, and Vj{vty is the group velocity of the state $* n (Co; E)(Qf a (( m+ i, E)) 

273 through the x-y plane at Co(Cm+i)- In the OBM scheme, group velocity is expressed as 



<!>Y(( ;Eyr L (<; ;E)$?(( ;E) 
$r°(Cm+i; E)^ R (( m+1 ; E)&r(c m+ i; E) 



(49) 



274 where L z is the length of the unit cell in the ^-direction. The proof of Eq. f f49|) is given in 

275 Appendix lAl 

276 The scattering wave function ^j(Ci] E) corresponding to the jth incident wave <£>* n (Co; E) 

277 is given by a linear combination of transmitted waves Qf a (Q;E) inside the right electrode 

278 (I > m + 1) with a transmission coefficient tij, i.e., 

N 

^■(0; E) = tij$T a ((i; E) = Q im (0; E) t 2j , ■ ■ ■ , t Nj f , (50) 



i=l 



16 



279 where Q tra (d;E) is an N- dimensional matrix that gathers right-propagating Bloch waves 

280 (lA(-E')l = 1 and Re(X(E)) > 0) and rightward-decreasing evanescent waves ([ A(-£7) | < 1) of 

281 ¥r((i;E),i.e., 



Q tra (Cf, E) = ¥r((f, E), $r°(0; E), ■ ■ ■ , S*F(C; E) 



From Eqs. (I40p and ( 150]) for I = m + 1, we have 



283 and then obtain 



(51) 



[t lh t 2j , • • • , t Nj f = iQ tra {C m+ u E)- l G r T {Cm+i, Co; £)r L (Co; E) (52) 



T = tQ tra (( m+ i; Ey^Cm+i, Co; E)T L (( ; £)Q m (Co; E). 



(53) 



284 Here, T is the transmission-coefficient matrix 



til ti2 
^21 ^22 



tlN 
tNN 



285 and Q m (Co; E) is the matrix defined by Eq. ( 1451) . Note that the expression 

Y^\U 3 \ 2 - =Tr[v- 1 T t V'T 



(54) 



(55) 



286 holds, with being a diagonal matrix whose elements are v\ ^Sy. By substituting Eq. ((53 
28? into the right-hand side of Eq. ( 1551) , we find 



Em 2 --' 1 - 



r L (Co; E)Q m (( ; i?)V- 1 Q m (Co; £) f r L (Co; E)G a T (Co, Cm+i; £) 



Q* TO (Cm+i; -E) f 1 V'(5* ra (Cm+i;-K) ^^(Cm+bCoj-E) 



(56) 



288 Here, 'Tr' stands for the trace, i.e., the sum of the diagonal matrix elements, and the cyclic 

289 property of the trace is used. From Eqs. (1151) . (jUJ), and (l5Tj) . the relations 



V = *W n (Co; i?) t r L (Co; £)Q m (Co; E) 

V = iL z Q tra (( m+l] Eyr R (( m+1 ; E)Q tra (( m+1 ; E) 



(57) 
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290 are derived, and then 



G 



h 




2e^ 
h 

h 



Tr ^(Co; E)Gt(( , Cm+ij E)T R (( m+1 ; E)G r T (( m+1 , ( ; E) 



Tr T L (( ; E)G r T (( , Cm+i! ^^(Cm+ii E)Gx(Cm+i, Co! -E) 



(58) 



29i are established 



Here, the advanced Green's function is 



Gr((k, Ci] E) = G r T ((i, ( k ; Ey . 



(59) 



292 Equation ( 1581) is a well-known formula^ in the NEGF formalism pioneered by Keldysh.- 

293 The equality of the last line in Eq. ( 158|) is also verified with a similar consideration of the 

294 case of incident waves incoming from the right electrode. One of the advantages of the 

295 Green's-function approach is that the conductance is calculated without the knowledge of 

296 well-defined asymptotic wave functions in the transition region. 

297 III. TRANSPORT PROPERTY OF BN RING CONNECTED TO CNT ELEC- 

298 TRODES 

299 To demonstrate the applicability of the RSFD NEGF method and the importance of the 

300 interpretation using scattering wave functions, the transport property of the C/BNNT where 

301 one carbon ring of (9,0) CNT is replaced with a BN ring is examined. Figure H] shows the 

302 computational model, in which the C/BNNT is sandwiched between the CNT electrodes. 

303 A valence electron-ion interaction is described using norm-conserving pseudopotentials 3 ^ 

304 generated by the scheme proposed by Troullier and Martins. 31 Exchange and correlation 

305 effects are treated within the local density approximation^ of the density functional theory. 

306 To determine the Kohn-Sham effective potential, we use a conventional supercell under a 

307 periodic boundary condition in all directions with a real-space grid spacing of ~ 0.24 A; the 

308 dimensions of the supercell are L x = 13.34 A, L y = 13.34 A, and L z = 4.32 A, where L x 

309 and L y are the lateral lengths of the supercell in the x- and y-directions perpendicular to 

310 the nanotube axis, respectively, and L z is the length in the ^-direction. Then we compute 
3n the scattering wave functions obtained non-self-consistently. It has been reported that this 
312 procedure is just as accurate in the linear response regime but significantly more efficient than 
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313 performing computations self-consistently on a scattering-wave basis.— The conductance is 

314 calculated using Eq. fl58|) . 

315 Figure [5] shows the conductance spectrum of the C/BNNT. To investigate states that 

316 actually contribute to electron transport, the eigenchannels are computed by diagonalizing 

317 the Hermit matrix T^T, where T is the transmission-coefficient matrix defined in Eq. ( 15^1) .— 

318 For reference, we plot in Fig. O the band structures of the CNT and C/BNNT in the periodic 

319 supercell used to obtain the Kohn-Sham effective potential. It is reported in the other DFT 

320 calculations that the band gap of the (9,0) CNT opens even though the (9,0) CNT is expected 

321 to be metallic by the zone folding of the tight-binding approximation;— the zero-conductance 

322 region around the Fermi level corresponds to the fundamental band gap of the (9,0) CNT. 

323 In addition, the bands of the C/BNNT indicated by (3 and 7 doubly degenerate as well as 

324 the bands of the CNT indicated by a at approximately the Fermi level. One can see in 

325 Fig. [5] that two channels actually contribute to the transport in the vicinity of the Fermi 

326 level and significant peaks due to the resonant tunneling through the dispersionless bands 

327 indicated by 7 are not observed in the conductance spectrum. Since a real-space picture of 

328 the wave function helps us to understand the transport phenomenon, the spatial behaviors 

329 of the C/BNNT and CNT states, which are relevant to the electron transport, are shown 

330 in Fig. [7J We plot the wave functions at V point because the symmetry in the x-y plane 

331 is insensitive with respect to the variation in k z in a one-dimensional Brillouin zone. For 

332 comparison, the behaviors of the wave functions indicated by a, (3, and 7 in Fig. [6] are also 

333 plotted in Fig. |HJ The wave function of the energetically dispersive bands of the C/BNNT 

334 (indicated by /3) shows threefold rotational symmetry with respect to tube axis, whereas 

335 that of the dispersionless bands (indicated by 7) shows fivefold symmetry. The energetically 

336 dispersive bands can contribute to electron transport because the spatial symmetry of the 

337 wave functions of the bands corresponds to that of the wave functions coming from the left 

338 CNT electrode and outgoing to the right CNT electrode indicated by a. In contrast, the 

339 symmetry of the dispersionless bands does not agree with those of the coming and outgoing 

340 waves, and thus the wave functions from the CNT electrodes are hardly connected to the 

341 wave function of the dispersionless bands, resulting in a small contribution to the electron 

342 transport. These results imply that the interpretation using the scattering wave function, 

343 which does not explicitly appear in the NEGF method, is important for the investigation of 

344 the transport properties of nanoscale systems. 
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34 5 IV. SUMMARY 



346 We have proposed an efficient procedure for RSFD NEGF calculations for obtaining 

347 the transport property of nanoscale systems. Since the number of grid points in real-space 

348 methods is larger than that of basis in the tight-binding approach, the computational cost to 

349 obtain the surface Green's functions and self-energy terms has been the bottleneck. Using the 

350 ratio matrices in the OBM method, the present procedure greatly reduces the computational 

351 cost to obtain the surface Green's functions and self-energy terms of the electrodes as well 

352 as the Green's functions of the whole system without loss of mathematical rigorousness. In 

353 addition, we proved that scattering wave functions, which provide us a real-space picture of 

354 the scattering process, can be obtainable by the present scheme. 

355 The transport property of the BN ring connected to the carbon nanotube electrode is 

356 investigated using the present method. By examining the rotational symmetry of wave 

357 functions at the matching plane, we found that states that are energetically dispersionless 

358 in a one-dimensional Brillouin zone hardly contribute to the electron transport because of 

359 the difference in the rotational symmetry of wave functions with respect to the tube axis. 

360 This result indicates that the real-space picture of the scattering wave function, which is 

361 not necessary to be taken into account to calculate the transport property in the NEGF 

362 formalism, helps us to interpret transport phenomena in the transition region. 

363 Since the OBM method is developed for the RSFD scheme, the present technique allows 

364 us to efficiently obtain quantities for the NEGF method in a real-space representation. In 

365 particular, the RSFD scheme of first-principles calculations is a method that has the advan- 

366 tage to scale with massively parallel architectures and has this potential without compromise 

367 on precision. Moreover, this scheme is free of problems concerning the completeness of the 

368 basis set such as in the methods using localized basis sets of either atomic orbitals or Gaus- 

369 sians. Therefore, the present procedure opens the possibility for executing large-scale RSFD 

370 NEGF calculations using massively parallel computers with a high degree of accuracy. 
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379 Appendix A: Group velocity 



380 For completeness, we introduce the expression for the group velocity in this appendix. 

381 Group velocity is written as 



dE 



Vn 



u\E)u{E), 



(Al) 



9 dk(E) 

382 where k(E) is a Bloch wave vector with E being the Kohn-Sham energy and u(E) is an 

383 N x x N y x iV^- dimensional columnar vector consisting of $j(^ M , E) of the Mth unit cell (see 

384 Fig. E]): 

u{E) = [^-(Cf; E), $,(C 2 M ; E), ■ ■ ■ , E)\\ (A2) 

385 The propagating wave obeys the Kohn-Sham equation 



H^ r (k(E),k\\)u(E) — Eu(E), 



(A3) 



386 where H^ r (k(E) is the N x x N y x A z -dimensional periodic Hamiltonian of the Mth unit 



cell, 



H™{k{E),k\ 



ff(Cf;fc||) B{C{ 



e -ik(E)L zB ^My 






(A4) 



388 Differentiating the eigenvalue equation of Eq. ( 1A3I) with respect to k(E) and multiplying 

389 u(Ey from the left-hand side of the resultant equation yields 



dH, 



M 



dE 



u{E)^u{E). 



(A5) 



From Eq. (IA4j) . one finds that 



dH. 



per 



dk(E) 








••• -tL z e- ik ^B((^ 
■■• 



iL z e lk W L *B((M) ■■■ 



(A6) 
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and 



-iL z e-^ L ^ 3 {^- E^BiO^jiC; E) 

-iL z ^+\E)^B{C)^AC-,E) 
-iL^C-^E^BiO^-.E) 
-iL^jiC^E^B^^iC^E). 

392 Here, the second step follows from the Bloch condition 

^(Ci M+1 ;^) = e^$ J (Cf;^), 



(A7) 



(A8) 



393 and the flux conservation is used in the last step. From Eqs. (IA1I) . (IA5I) and (IA7I) . group 

394 velocity is given by 



v g = iL z 



HCX- 1 ; eVb(C)^ 1 ; e) - $,(cf; E^BiO^AC- 1 ; e) 



(A9) 



395 where the normalization of the Bloch states is defined as ^ $j(£j M ; i£)t<&j(£j M ; E) — 1. 

396 Furthermore, making use of Eqs. fHTT) . ( |46|) . and (1471) . we have a simpler expression for 

397 group velocity in terms of the coupling matrix: 
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FIG. 1. Sketch of a system with the transition region intervening between the left and right semi- 
infinite crystalline electrodes. The dashed lines correspond to the borders of the partitioning of the 
Hamiltonian matrix in Eq. (PQ) and Fig. [2J whereas the dotted lines denote the borders of individual 
regions used in wave-function matching.— The case for A// = 2 and N z = 2m is illustrated as an 
example, where A// corresponds to the order of the finite-difference approximation for the kinetic 
energy operator. 
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FIG. 2. Partitioning of the Hamiltonian matrix H of Eq. ([T]), associated with the whole system 
sketched in Fig. [TJ Block-matrix elements Hi, B\ and Bw are the abbreviations of H(Q, fen), B(Q) 
and B(d,Q'), respectively. The partition lines are identical to those in Eq. ([T]). 
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FIG. 4. (color online) Computational model where one BN zigzag ring is sandwiched between (9,0) 
CNT electrodes. Large dark, small dark, and small light balls are N, C, and B atoms, respectively. 
The rectangle enclosed by broken lines represents the supercell used to evaluate the optimized 
atomic configuration and Kohn-Sham effective potential. 
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FIG. 5. Conductance spectra of C/BNNT as functions of energy of incident electrons. The solid 
curve represents conductance. Dashed, dotted, dash-dotted, and dashed double-dotted curves show 
channel transmissions. 



31 




r x r x 



FIG. 6. Electronic band structures of (a) CNT and (b) C/BNNT. The supercells include two and 
four (9,0) rings for CNT and C/BNNT, respectively. The zero of energy is the Fermi energy. 
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FIG. 7. Contour plots of wave functions, (a) and (b) Doubly degenerate bands indicated by (3 of 
C/BNNT in Fig. El (c) and (d) Doubly degenerate bands indicated by 7 of C/BNNT. (e) and (f) 
Doubly degenerate bands indicated by a of CNT. Negative values are indicated by dashed curves; 
the thick solid curve represents zero. The contour plot is separated by 7.41 xlO -4 electron/A 3 . 
The spheres represent the position of carbon atoms. 
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FIG. 8. Contour plots of doubly degenerate scattering wave functions at energy of Ep — 0.5 eV. 
Negative values are indicated by dashed curves; the thick solid curve represents zero. The contour 
plot is separated by 3.27 xlO -4 electron/ A 3 /eV. The spheres represent the position of carbon 
atoms. 
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